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Abstract 

O 

■ We propose a new Quantization algorithm for the approximation of 
inhomogencous random walks, which are the key terms for the valuation 
of CDO-tranches in latent factor models. This approach is based on a dual 
quantization operator which posses an intrinsic stationarity and therefore 
automatically leads to a second order error bound for the weak approxi- 

Qh ! mation. We illustrate the numerical performance of our methods in case of 

■ the approximation of the conditional tranche function of synthetic CDO 
£H , products and draw comparisons to the approximations achieved by the 

; j— | saddlepoint method and Stein's method. 

cr 
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i/-) ' 1 Introduction 

in ' 

\ In this paper we focus on the numerical approximation of inhomogeneous Ber- 

noulli random walks. 

Therefore, let (Q, F, (F), IP) be a filtered probability space on which we 
0^ \ define the inhomogeneous random walk 



X:=Y,aiZi, (1) 



i=l 

for some independent {0, l}-valued Bernoulli random variables Z% ~ B(j>i), Pi € 
(0, 1) and Qj > 0. 

The distribution of X plays a crucial role for the valuation of basket credit 
derivatives like CDO-tranches in latent factor models (see e.g. pQ or [5]). These 
are credit products, whose payoff is determined by the loss in large portfolios 
of defaultable credit under lyings. 
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Therefore assume that we have a portfolio of n defaultable credit names 
with notional amounts iVj and whose default times T{ are [Tt) stopping times, 
i = 1, . . . , n. Here, (J-t) stands for the observable filtration of the credit names. 
Moreover, we denote the fractional recovery of the i-th credit by Ri. 

Hence, the fractional loss of the portfolio up to time t is given by 

. ^ (l-R t )N m 

lt — Z^ n 1 ^<<}' (2) 

i=i 

where iV = X^=i ^Y? * s the total notional. 

Following the ideas of [6] and [3], the distributions of the default events 

< t} up to a fixed time t are driven under the risk-neutral probability 
measure by a common factor U (which we may assume w.l.o.g. as U([0, 1]) 
distributed) and some idiosyncratic noise £j. 

That means, that we assume that the events {r^ < t}, i = l,...,n are 
conditionally independent given a(U). 

Furthermore, we require the existence of a copula function F : [0, l] 2 — > 
[0, 1], such that p \— > F(p,u) is a non-decreasing, right continuous function for 
every u £ [0, 1] and 



f F(p,u)du = p, pe[0,l]. 

J 



Since it holds 



'{{ri < t}) = E(P({T; < t}\U)) = I F({Ti < t}\U = u) du 



J Q F(v({Ti<t}),u)du, 



we may interpret F^PUn < t}) ,u) as the conditional default probability 

P({n<t}\u = u). 

Typical choices for the function F are the standard Gaussian copula 



F(p, u) = $ 



3> -1 (p) - /9<l> _1 (li) 
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with common correlation parameter p, or the Clayton copula (cf. [S]). 

Thus, for a fixed time t, the risk-neutral conditional distributions of the 
portfolio losses It given the event {U = u} are driven by a random walk of type 
([1]) with oti := (1 — Rj)Ni/N and conditionally independent Bernoulli random 
variables Zi := l{ Ti <t} with parameters pi := F^P({Tj < t}) , uj. 

The cash flows of a (synthetic) CDO single tranche [a, b] with attachment 
points < a < b < 1 read as follows: 

The protections seller of the tranche [a, b] has to pay at each default time Tj 
which satisfies Z Ti € [a, 6] the notional of the defaulted name minus its recovery, 
i.e. 

(1- Ri)Ni. (default leg) 
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On the other hand he continuously receives a coupon payment of 

(premium leg) 



K Nl aM (t)dt, 



where k is the fair spread of the tranche, which is to be determined by arbitrage 
arguments. We denote by Nj a ' b \t) the outstanding notional of the tranche at 
time t, that is the notional amount of the tranche [a, b] which has not defaulted 
up to time t . 

Assuming a deterministic risk-free interest rate r and continuously com- 
pounding, we note that 

(1 - Ri)Ni 



N 



■t [aM {l Tl )=F^ b \n)-F^ b \r l -), 



where the tranche losses Fl a ' b ^ are defined as 
F\ aM (t) :=(l t - a)+ - (k - b) + = 



if l t < a 

It — a if a < It < b . 

b — a if It > b 



Hence, the discounted default payments accumulated up to maturity T maybe 
written as 

n n 

e~ rr > (i - Ri)Nil M (l n ) =NJ2 e~ rn [Fi aM in) - F\ aM (n-)] 1 T! < T 

i=l i=l 

= N [ T e~ rt Fl a > b] (dt). 
Jo 

Concerning the premium leg, the outstanding notional jvj a ' 6 ' (t) of the tranche 
[a, b] is given by 



Nj a ' b] (t) = N-[(b-a)-F^ b] (t)] 



N -(b-a) if l t < a, 
N ■ (b — l t ) if a < l t < 6, 
if l t > b 



so that the discounted coupon payments Ate rt N\ a,b \t) dt accumulate between 
and T to 



k-N I e~ rt [(b-a)-F\ aM (t)]dt. 



Under the risk-neutral probability measure both legs have to produce an 
equal present value, i.e. 



N 



f T e - rt Fj a ' b] (dt) = K ■ N f T e~ rt [(b -a)- F^ b] (t)] dt, 
Jo Jo 



so that taking (risk-neutral) expectation and processing an integration by parts 
yield the fair spread value k, namely 



[1 - e~ rT ] - Jo' e- rt EFj aM (t) dt 
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Here, the mathematical challenge consists in the computation of the expecta- 
tions EF^ a ' b \t). This leads, within the latent factor models, to the approxima- 
tion of the conditional expectations 

E(F} a ' b] (t)\U = u)= E((l t - a) + \U = u) - E((l t - b) + \U = u), 
since we have 

EFl aM {t)= / E(Fl aM {t)\U = u)du. (3) 
Jo 

As already announced, the conditional distribution of It is given by an inhomo- 
geneous random walk as defined in ([1]). 

We therefore focus in this paper on the approximation of the distribution 
of this type of random walks, the outer integral with respect to U in ([3]) can 
afterwards be approximated by standard quadrature formulae. 

For the usual applications n has a size of about 100, which is by far too large 
for an exact computation of the distribution of the random walk X, but still 
too small to get accurate approximations based on the asymptotics provided 
by limit theorems as n goes to oo. 

Moreover, we have to deal in this general setting with arbitrary coefficients 
ai, which destroy in general any recombining property of the random walk. As 
a consequence, no (recombining) tree approach can be implemented. 

So far, most approaches developed in the literature for the approximation 
of the conditional tranche expectation E(F^ a ' b ^\U) rely upon the saddle point 
method (cf. [7j) or an application of Stein's methods for both Gaussian and 
Poisson approximation (cf. [2]). 

Although based on completely different mathematical tools, both approaches 
suffer from the same lack of accuracy in the computation of 

n , 

E(J2<XiZi-K) 

i=l 

when the strike parameter K is "at-the-mean" , i.e. when Ym=i a iPi * s c l° se to 
K. From a theoretical point of view no control of the induced error is available. 
Finally, even if their numerical performances can be considered as satisfactory 
in most situations, these approximations methods are "static" : the "design" of 
the method cannot be modified to improve the accuracy if a higher complexity 
is allowed. 

The structure is as follows. In section 2 we introduce a new Dual Quantiza- 
tion scheme for the approximation of the inhomogeneous random walk ([1]). 
Moreover we establish error bounds for this approximation and discuss its 
asymptotic behaviour. Section 3 is devoted to the numerical implementation of 
this quantization scheme and its numerical performance. Finally, in section 4, 
we give a slight modification of this scheme to also capture the computation of 
sensitivities with respect to the probabilities pi and the coefficients ctj. 
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2 Approximation of inhomogeneous Random Walks 



We will focus in this section on the numerical approximation of the inhomoge- 
neous random walk 

n 

X = aiZi 

i=l 

for independent Z\ ~ B(pi), Pi £ (0, 1) and a>i > 0. 

An exact computation of the distribution of X is still not possible with 
nowadays computers, since in our cases of interest we have n ~ 100 and X has 
up to 2 n states. Hence we aim at constructing a random variable X with at 
most N <C 2 n states and which is close to X, e.g. K\X — X\ 2 is small. 

Due to the fact that there is no way to generate X directly, we have to 
construct approximations along the raondom walk 

X° = 0, 

X k = X k ~ l + a k Z k , k = l,...,n 

where the increment Z k is an ordinary Bernoulli random variable which is easy 
to handle. Clearly we have 

X = X n 

and of course this would work similarly in full generality, if X is a function of 
a Markov chain. 

Now suppose that we are equipped at each layer k with some grid T k = 
{x k , . . . , %% k } of size N k and a (possibly random) projection operator Ilr fe : 
R — ► Tfc, which maps the r.v.'s X k into T k . 

We then may state a recursive approximation scheme for X = X n as follows 

X° :=0 

X k ^nr^X^ + QfcZfc), k = l,...,n. 

This will be the main principle for constructing the approximation of X n . 
It remains to choose appropriate grids and projection operators IIr fe . Here, 
it will turn out that the obvious choice of nr fc as a nearest neighbor projection 
is not sufficient in this setting and we will have to develop a new approach. 

2.1 Quantization and Dual Quantization 

Regular Quantization In view of minimizing M\X — X\ 2 for a general r.v. 
X £ L 2 (P), the above problem directly leads to the well known quadratic 
quantization problem (cf. [3]) 

inf jE|X - X\ 2 : X r.v. with card{A > (^)} < n\ (4) 

at some level N £N. We will from now on call any discrete r.v. X Quantization 
and in particular if card{X(f2)} < N we call it N -Quantization. 
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In fact one easily shows that (U|) is equivalent to solving 



inf (E minlX - x\ 2 : T C R, cardir} < n\ , 
I zer > 

which means that Ilr would be chosen as a nearest neighbor "projection operator 
on r, i.e. 

where (C x (T)) x£ -p denotes a Borel-partition of R satisfying 
C X (T) C {^1: |e — ar| <min|£-y|}. 

Such a partition is called Voronoi- Partition (of IR related to V). 

In the one dimensional setting the Voronoi cell C Xi (T) generated by the 
ordered grid T = {x±, . . . ,xn} consists simply of the interval [ ^'^ , x 5 >+1 ] . 
Nevertheless we will use in this paper the more general notion of a Voronoi cell 
to emphasize the underlying geometrical structure and the fact that this can 
also be defined in a higher dimensional setting. 

One shows (see that the infimum in (J3J) actually holds as an minimum: 
there exists an optimal quantization X*' N (which takes exactly N values if X 
has infinite support). 

Concerning the approximation of an expectation, first note that for V := 
X(Q) we get 

EF(X) = ^2f(x)-F(X = x), (5) 

xer 

so that X in fact induces a cubature formula with weights ¥(X = x), x G T. 
This may provide a good approximation of KF(X), if X is close to the optimal 
solution of the quantization problem (jj]). 

For a Lipschitz functional F £ CLi p (R, R) we immediately derive the error 
bound 

\EF(X) - RF(X)\ < [F]up E\X - X\. 

If moreover F exhibits further smoothness properties, i.e. F £ C 1 (IR) with 
Lipschitz derivative, we may establish for a quantization X satisfying the sta- 
tionarity property 

E(X\X)=X, (6) 

a second order estimate (cf. [8]) 



\EF(X) - EF(X)\ < [F'] Lip E\X - X 
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Note that this stationarity property is always fulfilled if X is a solution to the 
optimal quantization problem 

In view of the Zador Theorem (Thm 6.2 in [1]), which describes the sharp 
asymptotics of the quantization problem @ as N goes to infinity, this leads to 
a quadratic error bound for an optimal quantization X* ,N of size N 

\EF(X) - EF(X* ,N ) | < C x ■ [F'] Lip ■ N~ 2 . 
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Unfortunately, in practice this stationarity property ([6]) is only satisfied if X 
is in some way optimized to "fit" the given distribution of X. This optimization 
is time-consuming and due to the complicated structure of X n not feasible in 
our case of interest. 

Hence we propose a (new) reverse interpolation operator to replace the 
nearest neighbor projection, which offers an intrinsic stationarity and therefore 
leads to a second order error bound without the need of adapting X to the 
exact distribution of X. 

Dual Quantization This alternative quantization approach for compactly 
supported random variables is based on the Delaunay representation of a grid 
r, which is the dual to its Voronoi diagram. Hence we will call this approach 
Dual Quantization. 

Suppose now to have an ordered grid T 



for a r.v. X with compact support included in [a, b] (Typically, [a, b] is the 
convex hull of the support of X). Moreover we introduce for convenience two 
auxiliary points xq := a and xn+i '■= b. 

The Delaunay tessellation induced by T then simply consists of the line 
segments Xj Xj+i, j = 0, . . . , N + 1, where we arbitrarily choose Xj Xj + \ to be 
the half-open intervals [xj,Xj+i) for j = 0, . . . ,N and xjy xjy+i as the closed 
interval [xn,xn+i]- This way we arrive at a true partition of the whole support 



To define a projection from X(Q) C [a, b] onto T, we will not just map any 
realization X(u>) to its nearest neighbor, but consider the two endpoints of the 
line segment Xj*Xj*+i into which it falls. 

We then perform a reverse random interpolation between these two points 
Xj* , Xj*+i in proportion to the "barycentric coordinate" 



i.e. we map X(uj) with probability A to xj* and with probability (1 — A) to 
Xj* + \ (see Figured]). 



a < x\ < X2 < ■ ■ ■ < xn < b 



of X. 



A : 



xj* +1 - X(u) 



Xj* +1 - Xj* 



Xj* 




A 



A 



X(u) 



Figure 1 



Reverse random Interpolation Operator J 



A formal definition of this operator is given as follows. 
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Definition 1. Let A ~ W([0.1]) be a r.v. on some probability space (fi, JF, P) 
and let T = (x±, . . . , xjv), xo := a, xn + i := b be an ordered of [a, b]. The Dual 
Quantization operator is denned by 



N 



Remark. Note that we can always enlarge the original probability space (0, JF, P) 
to ensure that A is defined on this space and is independent of any r.v. defined 
on the original space. Therefore we may assume w.l.o.g. that A is defined on 

(n,^,p). 

For £ € [a, 6] and Xj*Xj*+i denoting the line segment into which £ falls, we 

get 

njf(t) = Xj .)= Xr+1 ~ i and p^e) = x r+1 ) = i 



Xj*+l — Xj* Xj* + i — Xj* 

(7) 

so that J A satisfies the desired reverse interpolation property. 

As already announced, this Dual Quantization operator fulfills naturally a 
stationarity property: 

Proposition 1 (Stationarity). For any grid V = (xi, . . . , xjy) it holds 

E(J A (X)\X) = x. 

Proof. Let £ £ [a, b] and denote by Xj*Xj* + \ the line segment in which £ falls. 
Then note that 



N 

„ t(0 

„ * I ' T .■ i i — T, ; / I .T.' i i — T. ' I 



E(J r A (£) ) = E ^ (xil r , 7+1 -, \ (A) + x i+ i 1 r i (A) ) 1 



= xj. • P(Jr A (£) = xj.) + Zj*+i • P(Jr A (0 = ^*+i) 

= " ~ - + - Xj 

Xj* + i - Xj* V 

= £• 

The conclusion now follows from the independence of X and A which implies 
E(J r A (X)\X) = K(J r A (C))^ =x = X. 

□ 

Similar to the primal Quantization setting we then derive by means of the 
stationarity a second order estimate for the weak approximation of smoother 
integrands. 

Proposition 2. Let F G C 1 (M) with Lipschitz derivative. Then every grid F 
yields 

\EF(X) - EF{J A (X))\ < [F'] Lip ELY - jf(X)\ 2 . 
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Proof. Prom a Taylor expansion we derive 

\F(jf(X)) - F(X) - F'(X)(J T A (X) -X)\< [F'] Up \X - j£(X)\ 2 
so that the stationarity property (Proposition [T]) implies 

\E(F(j£(X))\X)-F(X)\ < [F'] Lip E(|X- j£{X)\ 2 \X). 
Taking expectations then yields the assertion. □ 

2.2 Application to the approximation of the inhomogeneous 
random walk 

2.2.1 The algorithm 

We are now in the position to design an approximation scheme based on Dual 
quantization in which the general projection operator r±r fc is replaced by the 
dual quantization operator J v k . Let IV . . . ,F n be some ordered grids. We set 

A k 

X° :=0 

(8) 

X k := J^{X k - x + a k Z k ), k = 1, . . . , n 

for Afe ~ U([0, 1]) i.i.d and independent of (^fc)o<fc<n- 

We wish to approximate EF(X n ) by its dually quantized counterpart EF(X n ) 

2.2.2 Error bound for the approximation of EF(X n ) 

Concerning the approximation power of the dual quantization scheme (JS|) for 
EF(X n ) with F € Cj-^ (R), we immediately derive from Proposition [2] the fol- 
lowing local error bound for any grid T k 

{EFiX*- 1 + a k Z k ) - EF(X k )\ < [F'] Lip E|(X fe-1 + a k Z k ) - X k \\ 

since X k = J^(X k ~ l + a k Z k ). 

As a matter of fact, the global error then consists of all the local insertion 
errors of the quantization operator along the random walk (X k )i< k < n . 

Theorem 1 (Global Error Bound). Let F £ C 1 (1R) with Lipschitz derivative. 
Then the Dual Quantization scheme ([SjJ related to the grids T k ,l < k < n, 
satisfies 

n 

\EF(X n )-EF(X n )\ < [F'] Lip ^^E\(X k ~ 1 + a k Z k ) — X k \ 2 

k=l 

n 

= [F'] Lip ^EliX^+akZk) - J^X^+akZkf. 

Proof. First note that it follows from Proposition [2] that for any o£l 
\EF(X + a) -EF(j£(X) + a)\ < [F'] Lip E\X - jf(X) | 2 . 
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Consequently, we get for any r.v. Z independent of X 



E 



F{X+Z)\Z 



E 



F(j£(X) + Z)\Z = z 



\EF(X+z) - EF{j£{X) + z)\< [F'] Up E\X - j£{X)\ 



and thus 

\EF(X + Z)- EF(j£{X) +Z)\< [F'] Up E\X - j£(X)f 
This finally yields 



l=k 



\EF(X n ) -EF(X n )\ <^EF(X k + ^ aiZi) — EF(X k ~ l + ^ a/Z/) 

k=l l=k+l 
n n 

= YyF(j^{X k ~ 1 +a k Z k )+ ^mZi) 

k=l l=k+l 

n 

- EF(X k - 1 +a k Z k + £ 

l=k+l 

n 

< [F']u P Y.m^+a^k) - J^{X k ~ l +a k Z h )\ 

k=l 

n 

= [F'] Up ^TElp^- 1 + a k Z k ) - X k f '. 



k=l 



□ 



2.3 Optimal choice of the grid V 

Let us temporarily come back to a static problem for an abstract random vari- 
able X with ¥(X £ [a, b]) = 1. In view of the second order estimate from 
Proposition [21 we arrive for a fixed number N G N at the optimization problem 

E\X-j£(X)\ 2 ^ inf . (9) 
rc[a,6],|r|<jv 

It is established in [11] that this infimum actually stands as a minimum. 
Hence optimal dual quantizers exists. Moreover the mean dual quantization 
error achieved by such an optimal grid differs from mean optimal quantization 
error of the primal quantization problem @ asymptotically only by a constant. 

Theorem 2 (Optimal rate (OH S])). Let X be a r.v. with V(X £ [a,b]) = 1 
and continuous density ip. Then it holds 



lim N 2 inf E\X-j£(X)\ 2 = 2 lim N 2 inf EminlX-d 2 = - 

Tc[a,6] AT— oo Tc[a,b] x^T 6 

|r|<7v |r|<7v 



W{z)f' 2 d 



Remark. This theorem about the asymptotics of the Dual Quantization problem 
can also be generalized to non compactly supported r.v.'s. and to non quadratic 
mean error (see |llj). 
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Given the formula of the gradient and the hessian of the optimization prob- 
lem @ with regard to T a Newton algorithm similar to the one described in [S] 
can be employed to construct numerically optimal dual quantization grids. 

Nevertheless, a straightforward alternative is to derive an (only asymptoti- 
cally optimal) dual grid from a grid which is optimal for the primal quantization 
problem Such grids are precomputed (cf. [TO]) and online available at 

www. quant izat ion. math-fi . com 

To transform these regular quantization grids into dual ones, we consider its 
midpoints, i.e. if yi, . . . , yjy denote an optimal grid for the primal quantization 
problem we simply define its dual grid 

x r .= y -i±mi, , i v i (io) 

This choice is motivated by the asymptotic formula of Theorem and its proof 
in [11], where exactly this midpoint rule establishes a connection between dual 
and regular quantization. Moreover, this connection allows to deduce the opti- 
mal rate for the dual quantizers from that for regular quantizers. 

Coming back to the problem of interest in this paper, the construction of 
optimal (primal or dual) grids for each X k , k = 1, . . . , n is clearly out of reach 
so that we have to make a "slightly" sub-optimal decision: we will choose grids 
which are optimal for a normal distribution matching the first two moments 
of X k , since such a A/"(//fc, <r|) distribution is close to X k for large values of k. 
Additionally we can restrict these grids to the convex hull of the support of X k , 

Moreover, our numerical observations even tend to confirm an optimal N~ 2 - 
rate for these sub-optimal grids. This emphazises again the importance of 
the intrinsic stationarity provided by the dual quantization operator in 
contrast to its primal counterpart, the nearest neighbor projection, where the 
stationarity only holds for grids specially optimized for the true underlying 
distribution, i.e. the r.v. X k in our case. 



3 Numerical implementation and results 
3.1 Numerical Implementation 

We now present numerical results and notes on the implementation of the Dual 
Quantization scheme ([8]) for the approximation of 



E(J2^Zi-K) , (11) 

i=i + 

by means of 

E(X n -K)+. (12) 

Concerning the second order estimate of Theorem [H the call function x i— > 
x+ clearly does not satisfy the assumptions of a continuously differentiable 
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function with Lipschitz derivative. Nevertheless we can replace x+ by <p e ( x ) '■= 
E(x + sY)+, where Y ~ A/"(0, 1) and e > to overcome this shortcoming. This 
function satisfies \tp e {x) — x+| < e, <^ £ G C°°(R) and < (p' E < 1. Furthermore 

_ * 2 

<p e writes <p e (x) = x ■ 3>(§) + "^^ e ^\ where is the distribution function of 
the standard normal distribution. 

We could imagine to compute E(A n — K) + using a backward dynamic pro- 
gramming formula based on ([8]). However such an approach is "payoff" depen- 
dent and consequently time-consuming since the computation needs to be done 
for many values of K as emphasized in the introduction. 

An alternative is to directly rely on the cubature formula 

N n +1 

E(X n ~ K) + = - K) + ■ P(X n = x?) 

i=0 

to approximate Here T n = {x™, . . . , x^v } is a dual grid of a normal 

distribution as described by (fT0|) in section 12.31 and we set Xq '. — 0, x^y^^ : — 

En 

The main task is then to compute the weights P(A n = x") for 1 < j < N n , 
which are given by the following forward recursive formula. 

Proposition 3. In the dual quantization scheme |2|) the weights P(X k = xf) 
satisfy 

N k +1 

P(A fc = xf ) = [(! " Pk) ■ Af (x^ 1 ) • P(A >fe ~ 1 = x k ~ u 



j=0 - 



+ Pk-\?(x*- 1 +a k )-F(X^ = x k J 



where 



if xf = x k * 



1 ~ x i-i( l 0> tfxf = x^ +1 
otherwise 



and j* := denoting the line segment, which satisfies £ € x^* x^„ +1 . 

Proof. We clearly have 

P(A fc = xf ) = P(X fe = x k \X k ~ l = x)" 1 ) • P(A >fc ^ 1 = x k ~ l ) 

3=0 
N k +1 

= ]T ¥{J^{X k ~ l + a fc Z fe ) = xf l^*" 1 = xj" 1 ) • P(A fc ^ = xj" 1 ). 
j'=o 

Since A^, are independent of A fc_1 , we derive 
P^pf*- 1 + a k Z k ) = xf (A^ 1 = x^- 1 ) = P(^ r A ;(x J fc ~ 1 + a k Z k ) = xf) 

= (l- Pk ).F(J^(x k - 1 )=x k ) 
+ p k -F(J^(x k - 1 + a k )=x k ) 
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so that finally ([7]) yields the assertion. 



□ 



Practical implementation From an implementational point of view we 
process the quantization scheme ([8]) and we start with 



X 1 



0. 



http://mathema.tician.de/dl/pub/pycuda-mit.pdf i.e. a grid Tq = {x®} = {0} 
and weight 

to? := ¥(X° = x?) = 1. 

To pass from time k — 1 to k, we suppose to have grids = {x\, . . . , x^} 

as described above. Additionally, we add the endpoints := and 2^+1 := 

Yli=i a i an d define := U {xq, x^r fc+1 }. Moreover we assume that the 
weights 



P(X fc_1 
have already been computed. 



fc-i 
w j 



), j = 0,...,JV fc _i + l 



We then could compute P(X fc = xf), 1 = 0,..., N k + 1 directly by means of 
Proposition However, this approach requires 2(iVfc + 2)(A^/ c _i + 2) evaluations 
of the barycentric coordinate Xf, which are not cheap operations, since each 
evaluation involves a nearest neighbor search to find the matching line segment 



Therefore it is more efficient to first iterate through the state space 

Tfc-i U (r fc _i + a k ) 



of the r.v. X + a k Z k . While computing for each £ £ r^-i U r^-i + a k its 
matching line segment x k - t x k - t+l in the grid T^, we directly update the weight 

vector (it!^) <i<jv fc +i a t positions Z = and I = + 1. 

This approach is given by Algorithm [1] and needs only 2(N k ~i + 2) nearest 
neighbor searches per layer k. 




Xj 1 + a k 



Figure 2: Weight-updating 
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Algorithm 1 Weight-Computation for the Dual Quantization scheme 



# Initialization 

r - {0} 



1 



10? <- 1 



for k = 1, . . . , n do 

,fc 



«,*«-(), j = 0,...,JV fc + l 

for J = 0, . . . , N k + 1 do 
# Case: Z fc = 



Find line segment Xj*Xj* +1 in which 1 falls 



Set 



fc A; 

J* 



^+=A-(l-p fc )-^- 1 
«>J.+i+= (l-A)-(l-p fc )-^ _1 

# Case: Z fe = 1 



Find line segment in which x k - 1 + falls 



Set 

<4+i+= (l-A)-pfe-^- 1 
end for 



end for 
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3.2 Speeding up the procedure 



3.2.1 Aggregation of insertion steps 

In view of the global error bound from Theorem [1] it is useful to reduce the 
number of grid insertion steps k • A natural way to do so, is to aggregate uq 
r.v. Zi into 

k-no 

z 'k = ^2 aiZi > k = l,...,n/n 

i=(k— l)no+l 

and then set 

X k = J^X"- 1 + Z' k ), k = l,...,n/n . 

E.g. with a choice of no = 2 we would insert a binomial r.v. with 4 states 
at every grid-point, but performing only 1/2 of the insertions. 

However, for this choice of no = 2 the overall number of nearest neighbor 
searches for the matching line segment Xj*xj* + \ remains the same as for no = 1. 



3.2.2 Romberg extrapolation 

An additional improvement of this method is based on the heuristic guess that 
the approximation of E_F(X) by EF(X n ) (see Proposition^ admits a higher 
order expansion 

EF(X) = EF(X n ) + kN~ 2 + o(iV~ 2 ). (13) 

We then may use quantization grids of two different sizes N\ <C N2 to cancel 
the second order term kN~ 2 in the above representation. 
This leads to the Romberg extrapolation formula 



Nf EF(A^) - JVf EF(X r j, 



EF(X) = ^ ^' J 2 + o(N^). (14) 



Although assumption (|13p is only of a heuristic nature, numerical results 
seem to confirm this conjecture (like for "regular" optimal quantization). 



3.3 Numerical experiments 

For the numerical results we implemented the above dual quantization scheme 
for grids of constant size 500 and 1000 in all layers k = 1, . . . , n. Regarding the 
Romberg extrapolation approach we applied the extrapolation formula (|14|) for 
sizes 100 and 500. 

As concerns methods to compare our approach to, we implemented a saddlepoint- 
point method (cf. [7]) and the Stein approach for a Poisson and Normal ap- 
proximation developed in [2]. 

We tested two typical situations: homogeneous and truly inhomogeneous 
Bernoulli random walks. 



15 




25 



30 



35 



40 



45 



50 



*Quantization@500 ♦Quantization@500-1 00 ■ 
■*• Stein Poisson —Stein Normal Saddle Point 



■Quantization@1000 



Figure 3: Absolute Errors for the call of various strikes (po = 0.05, a = 0.5). 



Homogeneous random walk. Let us start with a homogeneous Test- 
Scenario, i.e. all a.i are chosen equal to 1. Moreover we assume the pi to be 
a n sample of a log-normal distribution, which corresponds to the case of a 
Gaussian copula. Hence the parameters read as follows: 

• n = 100, 

• OLi = 1, 

• Pi=Po exp(<r£i - a 2 /2), & ~ jV(0, 1) i.i.d., 
with po G {0.05, 0.1, 0.2}, a = 0.5, 

Since this setting yields a recombining binomial tree, we can compute the 
exact reference values of 

E (i>-*0 + 

i=i 

for K G [0, 50] and plot the absolute errors as a function of the strike K to 
illustrate the numerical performances of the methods. This has been reported 
in Figures [3] to [SJ 

Inhomogeneous random walk I. To discuss a more realistic scenario, 
we present an inhomogeneous setting with on uniform distributed on the integers 
{1, 2, • • • , 10}, so that it is still possible to compute some reference values by 
means of a recombining binomial tree. The parameters read as follows 

• n = 100, 
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Figure 4: Absolute Errors for the call of various strikes (po = 0.10, a = 0.5). 




Figure 5: Absolute Errors for the call of various strikes (po = 0.20, a = 0.5). 



17 



0,035 



0,030 / \ 

0,025- /-- - V 

0,020 j \ 7 A 

0,015- I \ / \ 

0,010 / - \ / Y 

0,005- / \ / \ 

-o,oio / \ y 

0,015 : 

1 20 30 40 50 60 70 80 90 100 1 10 1 20 

*Quantization@500 ♦Quantization@500-1 00 -*-Quantization@1 000 
-►•Stein Normal —Saddle Point 

Figure 6: Absolute Errors for the call of various strikes (po = 0.05, a = 0.5). 
. on ~W{1,2,.-- ,10}, 

• Pi =p exp(cr& - a 2 /2), ~ AA(0, 1) i.i.d., 

€ {0.05,0.2}, a = 0.5, 

The numerical results are depicted in Figures and [71 Note that we have 
excluded the Stein-Poisson approach since this setting is already out of the 
Poisson-limit domain for po = 0.05 and consequently yield bad results. 

Inhomogeneous random walk II. Finally, we present a non-trivial case, 
where the a, are non-integer valued any more, i.e. we have chosen them to 
be U([0, 1]) distributed. Since in this setting the recombining property of a 
binomial tree is destroyed, we cannot compute the exact reference value any 
more. Therefore, we have chosen a grid of size iV = 10000 to compute a reference 
value, since such a large grid size yields in all former examples an absolute error 
less than 10~ 8 . To be more precise, the parameters has been chosen as follows: 

• n = 100, 

• Oi~W([0,l]), 

• pi = po exp(cr& - a 2 /2), & ~ J\f(0, 1) i.i.d., 
Po G {0.05,0.2}, a = 0.5. 

Since the Figures [8] and [9] are quite similar to those obtained in the first 
inhomogeneous setting (except a lower resolution), it seems very likely, that 
the former inhomogeneous setting is a very generic case to illustrate the general 
performance of the three tested methods. 



18 




Figure 7: Absolute Errors for the call of various strikes (po = 0.20, a = 0.5). 




Figure 8: Absolute Errors for the call of various strikes (po = 0.05, a = 0.5). 
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Figure 9: Absolute Errors for the call of various strikes (po = 0.20, a = 0.5). 



In all the above cases the quantization method remains very stable and 
outperforms even for a grid size of iV = 500 in nearly all cases the other tested 
methods. Only in the homogeneous setting and for very small probabilities pi, it 
cannot achieve the performance of the Stein-Poisson approximation. However, 
this excellence of the Stein-Poisson method in that particular setting is mainly 
caused by the fact, that the target distribution is an integer-valued one, as 
the Possion approximation is. Hence, these result are nontransferable to the 
inhomogeneous case. 

In the more complex inhomogeneous setting (Figures [6] and [7]) , we still ob- 
serve a strong domination of the quantization methods for small and moderate 
probabilities. Furthermore, in the case po = 0.2, we even get an error for the 
Romberg extrapolation with grid sizes 500 and 100, which is close to that of a 
1000-point quantization. 

Concerning the computational time for the processing of our Dual Quanti- 
zation algorithm, this approach is of course not as fast as the Stein's method, 
where one only needs to compute the two first moments of X and then evaluates 
the CDF-function of the standard normal distribution. To apply our scheme, 
we have to process at each layer k, < k < n at full grid similar to recom- 
bining tree methods. Nevertheless, the execution of Algorithm 1 implemented 
in C# on a Intel Xeon CPU@3GHz took for a grid size of N = 500 only a few 
milliseconds. Moreover, once the distribution of X n is established, we compute 
E(A n — K)+ for several strikes K (as needed in practical applications) in nearly 
no time. 

Finally, we want to emphasize, that this approach gives, through the free- 
dom to choose a larger grid size, a control on the acceptable error for the 
approximation . 
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4 Approximation of the Greeks 

Concerning the computation of sensitivities with respect to the parameters ai 
and pi, 1 < I < n, we consider / : E" x (0, l) n — > E, defined by 

(a,p) ^ /(a,p) := eQ^ a,Z 4 - . 

i=l 

We are now interested in the computation of J^- and J^-. 
Some elementary calculations reveal that for every Z € {0, . . . , n} 

%L = E(£ - (Jf - a,)) + - e(£ OiZi - K 
and 

so that our task consists of approximating the distribution of 

This can be achieved using a straightfoward adaption of the previous dual quan- 
tization tree, where we simply skip the l-th. layer. 
To be more precise, we set 



dtX := 
c\X~ k := I 



jJ^X*- 1 + a k Z k ) k + l, 
diX^ 1 k = l. 



Remark. This scheme can be processed simultaneously for all Z, 1 < / < n 
without increasing the number of nearest neighbor searches. 

Numerical experiments, which are not reproduced here, also confirm the 
good numerical performance of the dual quantization in this specific setting. 
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